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Abstract 

The  objective  of  this  project  was  to  develop  an  algorithm  to  accurately  determine 
atmospheric  density  using  simulated  GPS  data.  This  algorithm  is  designed  to  support  a 
future  USNA  Small  SateUite  mission.  Atmospheric  density  is  the  most  variable  factor  in 
orbit  propagation.  Thus,  the  uncertainty  in  density  generates  the  most  error  when 
predicting  a  satellite’s  future  position.  Numerous  models  have  been  developed  to 
account  for  the  variations,  but  more  accurate  models  are  needed. 

In  developing  the  algorithm,  SateUite  Tool  Kit®  (STK),  Analytical  Graphics  Inc.’s 
orbit  propagation  software,  was  used  to  generate  data  using  one  of  several  atmospheric 
models.  By  measuring  the  changes  in  the  satelUte’s  orbit  due  to  atmospheric  drag,  the 
density  was  accurately  calculated  to  within  1%  of  the  1976  Standard  Atmosphere  Model. 
To  vaUdate  the  algorithm,  the  density  output  was  compared  to  that  of  the  model  used  in 
STK®. 


The  USNA  SmaU  SateUite  Program  has  planned  to  design  and  place  a  sateUite  in 
low- Earth  orbit  (LEO)  with  a  GPS  receiver  on  board.  The  primary  mission  of  the 
sateUite  is  to  determine  density  in  the  upper  atmosphere.  Once  the  USNA  sateUite  is  on 
orbit,  the  algorithm  can  be  used  to  create  a  database  of  densities.  Other  smaU  satelUte 
programs  wiU  launch  similar  sateUites  to  generate  sufficient  data.  With  the  new 
atmospheric  density  data,  scientists  can  create  an  improved  atmospheric  model. 
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Introduction 

Density  variations  in  the  upper  atmosphere  are  common  and  current  density 
models  can  deviate  from  reahstic  conditions.  However,  with  sufficient  data,  scientists 
can  develop  newer,  more  accurate  models.  By  using  a  Global  Positioning  System  (GPS) 
receiver  on  board  a  sateUite,  one  can  obtain  data  in  the  form  of  time,  position,  and 
velocity.  With  these  data,  it  is  possible  to  calculate  the  instantaneous  atmospheric  density 
at  the  satelhte’s  position.  A  GPS  receiver  on  board  a  satellite  provides  the  capabihty  to 
acquire  position  and  velocity  data  continuously  and  inexpensively.  Thus,  density  can  be 
determined  throughout  the  satellite’s  orbit.  Current  methods  to  calculate  atmospheric 
density  use  terrestrial  observations  and  involve  averaging  through  one  or  more  orbits. 

This  project  investigates  the  process  in  calculating  atmospheric  density  from  GPS 
data.  A  program  is  used  to  process  the  position  and  velocity  data  to  compute  density  as  a 
function  of  altitude.  In  the  future,  this  program  can  be  used  to  process  GPS  data  on  board 
sateUites  specifically  designed  for  density  measurements.  With  sufficient  density  data 
from  numerous  sateUites,  future  scientists  can  use  these  data  to  develop  a  more  accurate 


model  of  atmospheric  density. 
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Background 

Atmospheric  density  models  are  extremely  valuable  in  astrodynamics 
applications.  The  major  use  of  density  models  is  for  orbit  propagation.  With  accurate 
models,  it  is  possible  to  determine  a  sateUite’s  future  position  with  a  high  degree  of 
certainty.  There  are  a  wide  variety  of  models  ranging  from  simphstic  static  ones  to 
complex  dynamic  models.  However,  all  of  these  models  have  substantial  errors  and 
deviate  from  actual  conditions  in  varying  degrees.  The  current  density  models  are  subject 
to  errors  of  approximately  15%.*  This  accuracy  is  insufficient  in  many  apphcations 
where  more  accurate  orbit  prediction  is  required.  For  example,  the  Air  Force  Space 
Command  aims  at  knowing  density  within  5%  of  the  actual  value. 

The  force  due  to  drag  is  one  of  the  major  perturbing  forces  on  satellites  in  low- 
Earth  orbit  where  atmospheric  density  is  the  most  variable  factor.  A  low- Earth  orbit 
(LEO)  sateUite  is  defined  as  a  sateUite  that  orbits  the  Earth  at  altitudes  ranging  from  300 
km  to  800  km  above  the  Earth’s  surface.  The  drag  force  causes  satelhtes  to  slowly  spiral 
towards  the  Earth.  Atmospheric  drag  makes  precise  orbit  predictions  difficult  because 
atmospheric  density  is  highly  variable. 

The  most  common  density  model  is  the  exponential  model.  Density  has  been 
shown  to  decrease  exponentially  with  height  in  accordance  with  the  following  equation  : 

p=Po*A  ^  )  (1) 

where  p  is  density,  Po  is  a  reference  density,  H  is  the  scale  height,  had  is  the  altitude  of  the 
sateUite  above  the  eUipsoidal  Earth  model,  and  ho  is  the  reference  altitude  corresponding 
to  the  reference  density.  Scale  height  is  a  function  of  the  atmosphere’s  composition  at  the 
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particular  altitude.  A  change  of  one  scale  height  results  in  a  decrease  of  density  by  a 
factor  of  - .  Although  there  is  no  standard  reference  for  the  exponential  model,  the 

e 

altitudes  and  densities  can  be  taken  from  existing  models  such  as  the  U.S.  Standard 
Atmosphere.  For  improved  accuracies,  the  different  reference  altitudes,  reference 
densities,  and  scale  heights  can  be  updated  frequently.  This  exponential  model  is  the 
most  basic  one,  but  is  not  adequate  for  high  precision  orbit  predictions. 

The  deviations  from  the  exponential  model  are  caused  by  three  main  influences. 
These  are  latitude  variations,  diurnal  variations,  and  the  1 1-year  solar  cycle.  Latitude 
variations  are  caused  by  the  Earth’s  oblateness.  Since  the  Earth’s  equatorial  diameter  is 
approximately  20  km  greater  than  the  polar  diameter,  a  satellite  in  a  circular,  inchned 
orbit  will  have  a  higher  altitude  at  higher  latitudes.  Therefore,  the  satellite  will  pass 
through  a  region  of  higher  density  closer  to  the  equator  than  at  higher  latitudes.  Diurnal 
variations  are  caused  by  the  Earth’s  rotation.  The  side  of  the  Earth  exposed  to  the  sun  is 
heated  and  the  atmosphere  expands.  At  a  given  altitude,  the  density  is  higher  because  the 
atmosphere  is  pushed  up  from  the  lower  altitudes.  On  the  nighttime  side,  the  atmosphere 
cools  and  contracts  and  the  density  is  less  than  the  daytime  side  at  a  given  altitude. 

EinaUy,  the  11-year  solar  cycle  affects  the  atmospheric  density  because  the  solar  radiation 
varies  with  time.  At  the  solar  maximum,  the  sun  releases  more  energy  and  the 
atmosphere  expands  as  the  Earth  receives  more  radiation. 

Atmospheric  models  have  been  created  to  account  for  these  variations.  Many 
atmospheric  models  in  the  past  have  been  constmcted  using  satellite  drag  analysis. 
However,  these  densities  were  derived  by  measuring  the  changes  in  the  orbital  elements 


over  at  least  one  orbit.  When  measuring  these  changes,  the  calculated  density  is  that  at 
perigee.  Yet  even  this  value  is  an  average  in  that  the  equations  are  averaged  over  one 
revolution.  The  resolution  of  the  location  where  the  density  was  calculated  is 
approximately  ±  20°  centered  at  perigee  due  to  the  averaging."^  Perigee  is  defined  as  a 
satellite’s  closest  point  of  approach  to  the  Earth,  while  apogee  is  the  farthest  point  in  a 
satellite’s  orbit,  as  seen  in  Figure  1. 


Figure  1:  Diagram  of  an  Elliptical  Orbit 

Another  problem  with  earher  atmospheric  models  is  the  sateUite  tracking 

capabihties.  According  to  Bminsma^,  in  the  past. 

Orbit  determination  did  not  allow  precise  estimation  of  the  sateUite 
position  at  a  given  time  essentiaUy  due  to  insufficient  tracking  capabUities 
as  weU  as  inaccurate  gravity  field  models.  Estimating  the  parameters  over 
an  orbit  of  several  days  was  the  only  way  to  reduce  random  errors  and 
geophysical  noise,  but  at  the  cost  of  lesser  resolution. 

By  using  GPS  receivers,  orbit  determination  can  be  extremely  accurate  as  well  as 
continuous.  Therefore,  it  is  possible  to  determine  the  atmospheric  density  at  a  given 
point  in  a  sateUite’ s  orbit  since  continuous  sateUite  position  and  velocity  measurements 
are  available.  Recent  improvements  in  accuracies  in  GPS  receivers  have  aUowed  for 
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many  new  applications  for  receivers.  There  are  many  types  of  GPS  receivers  available. 

One  GPS  receiver  that  is  well  suited  for  this  apphcation  is  the  Ashtech  G- 12  HDMA  GPS 
board.  The  position  accuracy  at  orbital  velocities  is  approximately  3.0  m  Circular  Error 
Probable  (CEP)  in  the  horizontal  direction  and  6.0  m  in  the  vertical  direction.  The 
Ashtech  receiver  has  a  velocity  accuracy  of  approximately  5  cm/s  at  orbital  velocities.^ 

This  accuracy  allows  for  precise  orbit  determination  and  GPS  receivers  can  be  used  to 
determine  atmospheric  density.  In  addition,  the  G- 12  HDMA  board  is  quahfied  to  orbital 
altitudes  and  velocities. 
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Theory 

A  satellite’s  orbit  is  defined  by  six  parameters  known  as  the  classical  orbital 
elements.  These  are  the  semimajor  axis  (a),  eccentricity  (e),  inclination  (i),  longitude  of 
ascending  node  (f2),  argument  of  periapsis,  and  tme  anomaly  (v ).  The  orbital  elements 
can  be  calculated  if  the  satelhte’s  position  and  velocity  are  known.  Figure  2  illustrates 
some  of  the  classical  orbital  elements,  while  the  terms  are  defined  in  Appendix  A.  In 
determining  density,  the  semimajor  axis,  eccentricity,  and  true  anomaly  must  be 
calculated. 


Figure  2:  Classical  Orbital  Elements^ 
V®  :  True  Anomaly 
®  :  Argument  of  Periapsis 
£1  :  Longitude  of  Ascending  Node 
e  :  Eccentricity  Vector 
i :  Inclination 
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The  semimajor  axis  defines  the  size  of  the  orbit.  Given  position  and  velocity,  the 
semimajor  axis  can  easily  be  calculated.  The  specific  mechanical  energy  of  the  orbit  is 
defined  in  equation  (2)  by  the  sum  of  the  specific  potential  and  kinetic  energies  of  the 
sateUite: 


where  ^  is  the  specific  energy,  v  is  the  magnitude  of  the  velocity  vector,  r  is  the 
magnitude  of  the  position  vector,  and  p.  is  the  Earth’s  gravitational  parameter.  The 
Earth’s  gravitational  parameter  is  defined  as  the  product  of  the  Earth’s  mass  and  the 
gravitational  constant  in  Newton’s  Universal  Eaw  of  Gravitation.  As  shown  by 

Q 

VaUado  ,  once  the  energy  is  calculated,  the  semimajor  axis  {a)  can  be  determined  by 
using  the  following  relationship: 


(3) 


The  eccentricity  measures  the  shape  of  the  orbit.  An  orbit  with  an  eccentricity  of 
0  defines  a  circular  orbit  and  as  the  eccentricity  increases  towards  1,  the  orbit  becomes 
more  elliptical.  The  eccentricity  vector  (e  )  is  defined  as: 


The  magnitude  of  the  eccentricity  vector  represents  the  shape  of  the  orbit,  while  the 


vector  points  in  the  direction  of  peiiapsis. 
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The  final  orbital  element  necessary  in  density  calculations  is  the  tme  anomaly. 

This  is  an  angle  measured  from  the  eccentricity  vector  to  the  position  vector  as  defined 
by  the  following  equation: 


(5) 


where  v  is  the  tme  anomaly,  e  is  the  eccentricity  vector  and  f  is  the  position  vector. 

Atmospheric  drag  is  a  function  of  the  atmospheric  density,  the  satellite’s 
coefficient  of  drag,  cross  sectional  area,  mass,  and  velocity  relative  to  the  atmosphere. 
The  acceleration  due  to  drag  is  given  by  the  equation: 


1  Cn  *A  2  ^  ] 

Q  — _ Q. _ * _ — 

“drag  o  K  ''  re/  I-  I 

2  m  V 


In  equation  (6),  Cd  represents  the  coefficient  of  drag,  A  is  the  cross-sectional  area 


of  the  sateUite,  p  is  the  atmospheric  density,  and  Vrei  is  the  sateUite’s  velocity  relative  to 


the  Earth’s  atmosphere.  The  sateUite’s  baUistic  coefficient  is  defined  as - .  As  a 

C^-*A 

sateUite  in  LEO  orbits  the  Earth,  it  experiences  a  force  from  atmospheric  drag  in  the 
direction  opposite  the  velocity  vector. 

By  measuring  the  change  in  the  semimajor  axis,  it  is  possible  to  calculate  the 
density.  In  an  equation  developed  from  the  Gaussian  variation  of  parameters,  VaUado^ 
derives  an  equation  for  the  time  rate  of  change  of  the  semimajor  axis: 
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In  this  equation,  the  time  rate  of  change  of  the  semimajor  axis,  — ,  is  a  function  of  the 

dt 

baUistic  coefficient,  density,  relative  velocity,  eccentricity,  tme  anomaly,  and  mean 
motion.  A  derivation  of  equation  (7)  can  be  found  in  Appendix  D.  Mean  motion,  a 
function  of  the  semimajor  axis  and  the  Earth’ s  gravitational  parameter,  is  defined  as: 


Solving  equation  (7)  for  density. 


The  terms  in  equation  (9)  can  be  calculated  from  position  and  velocity  information  with 
the  exception  of  the  baUistic  coefficient  and  the  time  rate  of  change  of  the  semimajor 
axis.  The  baUistic  coefficient  can  be  variable  depending  on  the  mass  of  the  sateUite  and 
its  orientation.  If  the  sateUite  has  a  propulsion  system,  the  mass  wiU  decrease  as  fuel  is 
expended,  thus  changing  the  baUistic  coefficient  with  time.  The  cross-sectional  area  can 
change  if  the  spacecraft’s  orientation  is  changing  with  respect  to  the  atmosphere  and  the 
direction  of  motion.  To  eliminate  these  variables,  the  sateUite  used  to  determine  the 
density  must  have  a  constant  mass  as  well  as  a  constant  cross  sectional  area.  Therefore, 
the  sateUite  would  be  designed  with  no  propulsion  system.  To  eliminate  orientation 
problems,  the  sateUite  would  be  designed  to  be  spherically  symmetric.  FinaUy,  the 
coefficient  of  drag  can  be  approximated  as  2.2  for  spherical  sateUites.^'^  The  time  rate  of 
change  of  the  semimajor  axis  is  the  most  important  variable  in  equation  (9).  To  obtain 
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the  term  — ,  the  numerical  derivative  of  the  semimajor  axis  with  respect  to  time  must  be 
calculated. 

Several  assumptions  were  made  when  calculating  the  atmospheric  density.  In 
equation  (9),  the  time  rate  of  change  of  the  semimajor  axis  must  be  isolated  due  to  drag 
only.  However,  the  other  major  perturbing  force  on  a  sateUite  in  LEO  is  due  to  the 
Earth’s  geopotential.  The  Earth  is  shaped  as  an  oblate  spheroid  thus  its  gravitational  field 
is  not  a  perfect  sphere.  The  Earth’s  asymmetric  gravitational  field  causes  perturbations  in 
a  sateUite’s  orbit  and  it  affects  the  semimajor  axis  with  periodic  variations.  To  determine 
the  change  in  the  semimajor  axis  due  to  drag,  I  assumed  the  effects  on  the  semimajor  axis 
due  to  the  Earth’s  geopotential  and  drag  were  independent.  Thus,  to  determine  the 
changes  in  the  semimajor  axis  due  to  drag,  the  effects  of  the  geopotential  could  be 
algebraically  subtracted  from  the  data.  This  assumption  was  found  to  be  incorrect  as  seen 
in  the  Results  and  Discussion  section;  however,  by  limiting  the  length  of  time  in  which 
the  two  were  subtracted,  the  differences  were  reasonable.  Thus,  the  assumption  was 
adequate  for  this  apphcation. 

Another  assumption  is  that  the  only  major  forces  on  the  sateUite  are  drag  and  the 
geopotential.  Outside  forces  such  as  third- body  effects  and  solar  radiation  are  assumed  to 
be  much  smaUer  than  these  forces  and  do  not  affect  the  sateUite.  This  is  reasonable  in 
that  in  low- Earth  orbit,  in  this  case,  at  altitudes  below  500  km,  these  forces  are  small. 

This  can  be  seen  in  Eigure  3  for  an  orbit  with  a  semimajor  axis  of  6878  km. 

The  third- body  effects  of  the  sun  and  moon  as  weU  as  the  effect  of  solar  radiation 
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Figure  3:  Periodic  Effects  on  the  Semimajor  Axis  of  Third-Body  Perturbations 
and  Solar  Radiation  Compared  with  Atmospheric  Drag 

are  pericxiic  and  much  smaller  than  the  effect  of  atmospheric  drag.  Therefore,  the 

perturbations  which  affect  the  semimajor  axis  the  most  are  atmospheric  drag  and  the 

Earth’s  asymmetric  gravity  field.  The  effects  of  the  Earth’s  gravity  field  will  be  shown  in 

the  Results  and  Discussion  section  of  this  report. 

The  final  assumption  was  that  the  atmosphere  was  rotating  with  the  Earth.  This 
assumption  is  reasonable  in  that  the  U.S.  Standard  Atmosphere  models  the  atmosphere  to 
rotate  with  the  Earth.  In  addition,  the  atmosphere  is  not  stationary  and  rotates  at 
velocities  up  to  the  Earth’s  rotational  velcxity.  Equations  developed  by  King-Hele^  *  were 
used  to  model  the  rotating  atmosphere  in  my  algorithm.  These  equations  can  be  found  in 


Appendix  E. 
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Methodology 

I  developed  several  steps  in  my  algorithm  to  calculate  atmospheric  density.  These 
steps  were: 

1 .  Rotate  the  GPS  data  into  the  Geocentric  Equatorial  Frame 

2.  Calculate  orbital  elements  from  GPS  data 

3.  Calculate  the  sateUite’s  relative  velocity. 

4.  Isolate  the  effects  of  drag  on  the  semimajor  axis 

5.  Calculate  the  time  rate  of  change  of  the  semimajor  axis 

6.  Calculate  the  atmospheric  density 

First,  GPS  data  must  be  processed  to  determine  the  classical  orbital  elements.  I 
used  Satelhte  Tool  Kit®^^  to  simulate  the  GPS  data  in  the  form  of  position  and  velocity  in 
the  Earth- Centered,  Earth- Fixed  coordinate  frame.  In  addition,  time  and  altitude  are  also 
included  in  the  simulated  GPS  data.  These  data  included  effects  of  the  Earth’s 
geopotential  and  atmospheric  drag.  The  disturbing  potential  for  the  Earth’s  nonspheiical 
shape  is  given  by  the  following  equation: 


Pim  [sin((|)^^7-)]{C/„ 


Im 


(10) 


where  q  is  the  Earth’s  gravitational  parameter,  r  is  the  magnitude  of  the  sateUite’s 
position  vector,  is  the  Earth’s  radius,  P,^[sin((l)^^^)]  is  the  associated  Eegendre 
function,  is  the  sateUite’s  latitude,  is  the  sateUite’s  longitude,  and  and 
Si^  are  coefficients  that  model  the  Earth’s  spherical  harmonics  .  In  equation  (10),  l  and 
m  describe  the  degree  and  order  of  the  geopotential  model.  Satellite  Tool  Kit®  uses  the 
Joint  Gravity  Model  (JGM)  2  as  its  geopotential  model.  In  my  analysis,  I  used  a  70x70 
(/  X  m )  geopotential  model  and  the  1976  Standard  Atmosphere  to  generate  the  GPS  data. 
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Initially  a  time  step  of  60  seconds  was  used  over  one  orbit,  yielding  approximately  90 
data  points. 

With  the  simulated  GPS  data,  position  and  velocity  were  rotated  from  the  Earth- 
Centered,  Earth- Fixed  frame  to  the  inertial  Geocentric  Equatorial  frame.  In  addition,  the 
velocity  of  the  satelhte  relative  to  the  atmosphere  was  calculated.  Following  this,  the 
classical  orbital  elements  were  calculated  at  each  time  step.  These  data  show  how  the 
orbital  elements,  mainly  the  semimajor  axis,  vary  with  time.  However,  these  data  include 
effects  of  the  geopotential  with  drag.  For  density  analysis,  the  changes  in  semimajor  axis 
with  respect  to  time  due  to  drag  are  needed. 

To  obtain  the  changes  in  the  semimajor  axis  due  to  drag  only,  I  ran  two  more 
STK®  simulations.  The  first  was  the  same  orbit,  however,  only  due  to  the  geopotential 
effects,  while  the  second  was  the  two-body  case.  The  two-body  case  is  defined  by 
assuming  the  Earth  and  the  satellite  are  point  masses.  The  other  assumption  in  the  two 
body  case  is  that  the  satellite  is  not  subjected  to  any  forces  other  than  gravity  of  the  Earth. 

With  each  of  these  two  runs,  STK®  output  the  semimajor  axis.  I  then  subtracted  out  the 
semimajor  axis  effects  due  to  the  geopotential  from  the  semimajor  axis  due  to  the 
geopotential  and  drag.  Finally,  to  obtain  the  semimajor  axis  due  to  drag  only,  I 
algebraically  added  the  semimajor  axis  of  the  two-body  case. 

The  next  step  in  the  process  is  to  determine  the  time  rate  of  change  of  the 
semimajor  axis.  To  do  so,  I  used  the  MATEAB®*"^  function  “polyfit”  to  model  the 
variation  of  the  semimajor  axis  with  respect  to  time  as  a  polynomial.  With  this 
polynomial,  I  took  the  numerical  derivative  to  obtain  the  time  rate  of  change  of  the 
semimajor  axis.  With  this  information,  aU  of  the  values  in  equation  (9)  are  known. 
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Density  as  a  function  of  time  and  more  importantly,  density  as  a  function  of  altitude  can 
be  determined. 
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Results  and  Discussion 

In  conducting  my  research,  I  was  able  to  calculate  atmospheric  density,  obtaining 
percent  errors  consistently  below  1%  when  compared  to  the  1976  Standard  Atmosphere 
model.  However,  in  the  process  there  were  several  phenomena  I  had  not  anticipated. 
These  were  the  geopotential  effects,  latitude  effects,  and  the  atmospheric  rotation.  Once 
my  algorithm  was  generating  density  values  corresponding  to  the  input  model,  I  added 
random  noise  into  the  GPS  data.  I  also  included  complex  density  models  such  as  the 
Harris -Priester,  Jacchia- Roberts,  Jacchia  1960,  and  Jacchia  1971  models  in  the  STK® 
simulations.  With  these  models,  diurnal  variations,  geomagnetic  variations,  and  the 
effects  of  the  1 1 -year  solar  cycle,  were  examined  with  respect  to  changes  in  the 
atmospheric  density. 

I  produced  plots  of  density  vs.  altitude  for  several  different  orbits.  These  included 
circular  equatorial,  circular  inchned,  eccentric  equatorial,  and  eccentric  inclined  orbits. 

In  addition,  the  semimajor  axis  was  varied  to  account  for  altitudes  between  250  km  and 
400  km.  Figure  4  shows  a  plot  of  density  vs.  altitude  for  the  circular  equatorial  case,  with 
an  initial  semimajor  axis  of  6678  km.  In  the  plots  of  density  vs.  altitude,  the  “Drag 
Analysis”  case  represents  the  output  of  my  algorithm.  Figure  5  plots  the  percent  error  vs. 
time  for  the  same  case.  The  errors  with  the  circular  equatorial  case  are  under  1%  and 
remain  below  the  1%  bound  for  one  orbit.  Corresponding  plots  for  inchned  orbits  are 
presented  in  Appendix  B.  The  percent  error  in  the  calculated  density  for  an  inclined  orbit 
began  within  the  1%  error  bound,  but  the  inchned  orbits  were  subject  to  an  increasing 
error  in  the  calculated  density  as  time  progressed.  The  increasing  errors  were  due  to  the 
assumption  that  effects  of  the  geopotential  on  the  semimajor  axis  were  independent  of  the 
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effects  of  drag.  The  variations  in  the  semimajor  axis  due  to  the  geopotential  are  a 
function  of  inclination  because  of  the  shape  of  the  Earth.  The  following  plots  (Figures  4 
and  5),  along  with  others  found  in  Appendix  B,  represent  the  final  output  of  my 
algorithm. 


Figure  4:  Calculated  Density  and  Reference  Density  vs.  Altitude: 
Circular  Equatorial  Orbit 


Geopotential  Effects 

The  first  major  finding  was  that  changes  in  the  semimajor  axis  due  to  drag  and  the 
geopotential  were  not  independent  and  that  the  initial  assumption  was  incorrect.  Analysis 
using  STK®  demonstrated  that  there  was  a  couphng  between  the  two  forces  as  shown  by 
Figure  6.  This  plot  shows  the  semimajor  axis  of  a  sateUite  in  a  circular  orbit  with  an 
inclination  of  30°.  The  effects  of  the  geopotential  are  periodic  and  cause  the  semimajor 
axis  to  vary  in  magnitude.  The  straight  line  in  Figure  6  plots  the  semimajor  axis  in  the 
two-body  case,  where  the  Earth  is  modeled  as  a  perfect  sphere,  and  thus,  has  a  uniform 
gravity  field.  The  sinusoidal  tine  shows  the  periodic  effects  of  the  geopotential  on  the 
semimajor  axis. 


Semimajor  Axis  vs.  Time 


Figure  6:  Periodic  Variations  in  the  Semimajor  Axis  Due  to  the  Earth’ 

Geopotential 
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The  coupling  of  the  changes  in  semimajor  axis  due  to  the  geopotential  and  drag 
can  be  seen  in  Figures  7  and  8.  Figure  7  plots  the  changes  in  the  semimajor  axis  due  to 
the  geopotential  as  weU  as  due  to  the  geopotential  combined  with  drag.  The  plot  of  the 
changing  semimajor  axis  due  to  the  geopotential  and  drag  combined  shows  that  there  is  a 
decrease  as  time  progresses. 


Semimajor  Axis  vs.  Time 


Figure  7:  Variation  in  the  Semimajor  Axis  due  to  the  Earth’s  Geopotential 
and  Geopotential  and  Drag 

The  coupling  can  be  seen  in  Figure  8,  where  the  semimajor  axis  that  was  calculated 
decreases,  but  has  periodic  variations  as  it  decreases.  The  periodic  effects  are  due  to  the 
coupling  with  the  geopotential.  As  the  sateUite’s  semimajor  axis  decreases,  the  density 
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increases,  and  therefore  the  rate  at  which  the  semimajor  axis  decays  due  to  drag  increases 
as  weU.  Therefore,  the  rate  at  which  semimajor  axis  changes  due  to  drag  is  not  constant. 

As  time  progresses,  the  differences  between  the  pure  drag  case  and  the  approximated 
drag  case  increase.  The  pure  drag  case  is  where  the  only  perturbation  input  into  STK® 
was  atmospheric  drag.  In  the  approximated  drag  case,  the  semimajor  axis  was  calculated 
by  algebraically  subtracting  the  geopotential  only  case  from  the  case  of  the  geopotential 
combined  with  drag. 


Semimajor  Axis  vs.  Time 


Figure  8:  Variation  in  the  Semimajor  Axis  Due  to  Drag:  Calculated 

and  Pure  Drag 

To  correct  for  the  error  growth,  I  determined  that  by  limiting  the  length  of  time 
for  the  data  collection  to  half  of  a  period,  the  differences  between  the  two  cases  are 
within  reasonable  bounds  of  1%.  Once  the  satelhte  has  made  one  half  of  a  revolution,  the 
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new  semimajor  axis  is  input  into  the  STK®  simulation  so  that  the  deviations  do  not 
increase  to  an  unreasonable  magnitude.  The  couphng  between  the  effects  of  the 
geopotential  and  drag  are  responsible  for  the  increasing  errors  in  my  density  calculations. 
In  the  circular  equatorial  case,  the  low  percent  error  of  1%  can  be  attributed  to  the 
neghgible  geopotential  effects.  Since  the  inclination  is  zero,  the  magnitude  of  the 
variations  in  semimajor  axis  are  insignificant  and  thus  it  is  not  a  major  source  of  error. 

For  inchned  orbits,  the  errors  increase  with  time  as  shown  in  Figure  9.  This  plot  is  of  a 
circular  orbit  with  an  inclination  of  30°.  Goals  for  the  errors  in  the  density  calculation 
were  set  to  be  less  than  or  equal  to  1%.  Table  1  indicates  the  values  of  the  errors  with 
time.  Errors  remain  at  a  reasonable  level  of  1%  through  approximately  half  a  revolution. 
Past  this  point,  the  errors  continue  to  grow.  Past  one  revolution,  the  errors  become 
inadequate  for  accurate  density  analysis.  Therefore,  the  orbital  elements  must  be  updated 
at  a  minimum  of  half  a  revolution  to  ensure  errors  of  1%  or  less. 


Table  1:  Long-Term  Errors  in  Density  Calculations 


Time 

Error  (%) 

1/4  Orbit 

<  0.5 

1/2  Orbit 

1 

3/4  Orbit 

3 

1  Orbit 

5 

2  Orbits 

19 
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Percert  Enrtx  vs.  Time 


Figure  9:  Long  Term  Percent  Error:  Circular  Inclined  Orbit  (Inclination  30°) 
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Latitude  Effects 

Latitude  effects  on  density  calculations  played  a  large  part  in  the  accurate  density 
calculations.  Initially,  to  determine  the  sateUite’s  altitude  I  subtracted  the  magnitude  of 
the  satellite’s  position  vector  from  the  Earth’s  radius.  Thus,  the  calculated  altitude 
represented  that  for  a  spherical  Earth  rather  than  the  actual  oblate  Earth.  As  seen  in  a 
circular  inclined  case  in  Eigure  10,  the  density  corresponds  near  perigee  and  apogee  (300 
and  283  km),  where  the  satellite  is  at  the  equator.  During  the  course  of  an  orbit,  the 
satellite  begins  at  an  apogee  of  300  km,  and  decreases  in  altitude  to  approximately  283 
km  at  perigee.  The  sateUite  then  returns  to  its  initial  altitude  of  300  km.  In  this  case,  the 
satellite’s  inclination  is  30°.  However,  as  the  sateUite  increases  in  latitude,  the  actual 
altitude  is  higher  than  the  calculated  altitude.  Therefore,  the  calculated  density  is  lower 
between  perigee  and  apogee  than  the  actual  values. 

By  using  STK®  generated  altitudes  that  took  the  Earth’s  oblateness  into  account, 
the  calculated  density  corresponded  with  the  U.S.  Standard  Atmosphere  model.  The 
addition  of  the  altitude  information  is  reasonable  because  GPS  has  the  capabiUty  of 
outputting  the  altitude  of  the  receiver,  where  the  altitude  is  referenced  to  an  accurate 
ellipsoidal  Earth  model. 


Density  (kgftn  ) 
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vio  "  Density  vs.  Altitude 


Figure  10:  Latitude  Effects  on  the  Calculated  Density:  Circular  Inclined  Orbit 
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Atmospheric  Rotation 

The  rotating  atmosphere  also  affects  density  calculations.  Although  the  value  of 
the  atmosphere’s  velocity  is  small  when  compared  to  the  satellite’s  velocity,  the  effect  on 
the  calculated  density  can  cause  an  error  as  large  as  12%.  The  effects  of  the  atmospheric 
rotation  can  be  seen  in  the  following  three  plots,  Figures  11,  12,  and  13.  These  three 
density  calculations  assume  a  non- rotating  atmosphere.  Figure  1 1  shows  density 
calculations  for  a  direct  orbit,  where  the  satellite  has  an  inclination  of  0°  .  The  plot  shows 
a  consistent  error  of  approximately  12%  below  the  reference  densities.  In  the  case  of  a 
direct  orbit,  the  sateUite’s  relative  velocity  is  lower  than  its  inertial  velocity.  Thus,  the 
calculated  densities  are  lower  than  the  reference  densities. 


Figure  11:  Calculated  Density  and  Reference  Density  vs.  Altitude: 

Circular  Equatorial  Orbit 
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The  next  plot,  Figure  12,  shows  density  calculations  from  a  polar  orbit,  which  has 
an  inclination  of  90°.  In  this  mn,  the  calculated  density  corresponds  with  the  U.S. 
Standard  Atmosphere  densities,  however  the  errors  increase  as  time  progresses.  It  is 
assumed  in  this  case  that  the  sateUite’s  velocity  vector  is  perpendicular  to  the 
atmosphere’s  velocity,  and  does  not  have  an  effect  of  the  density  calculation. 


Density  vs.  AMude 


Figure  12:  Calculated  Density  and  Reference  Density  vs.  Altitude: 

Circular  Polar  Orbit 

Finally,  Figure  13  shows  a  retrograde  orbit  with  an  inclination  of  180°.  In  this 
case,  the  sateUite  is  going  “against”  the  wind  and  the  relative  velocity  is  lower  than  its 
inertial  velocity.  This  plot  has  an  error  of  12%  when  compared  to  the  reference  densities; 


however  unlike  the  direct  orbit,  the  calculated  densities  are  higher  than  the  reference 
densities. 
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Y  1 Q '  ’  Density  vs.  Altitude 


Figure  13:  Calculated  Density  and  Reference  Density  vs.  Altitude: 

Circular  Equatorial  (Retrograde) 

These  three  plots  demonstrate  that  there  is  a  significant  effect  on  the  density 
calculations  and  that  is  a  function  of  the  sateUite’s  inclination.  Since  the  direct  orbit  and 
the  retrograde  orbit  produce  errors  of  the  same  magnitude,  but  different  signs,  the 
atmosphere’s  velocity  influences  the  calculations  considerably.  Once  the  atmosphere’s 
velocity  was  accounted  for,  the  errors  were  significantly  reduced. 
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Complex  Atmospheric  Models 

There  are  many  atmospheric  models  that  include  complex  features  of  the 
atmosphere  such  as  the  effects  of  diurnal  variations,  varying  solar  flux,  and  geomagnetic 
activity.  The  algorithm  developed  illustrates  these  variations  quantitatively.  Although 
the  calculated  density  values  are  not  vahdated  with  the  corresponding  models,  the 
patterns  are  consistent  with  expected  trends.  Atmospheric  density  reaches  its  maximum 
at  approximately  1400  local  solar  time^^.  This  trend  is  illustrated  in  Figure  14.  This 
figure  plots  four  mns  where  a  different  atmospheric  model  was  used  in  the  orbit 
propagation.  The  orbit  used  in  Figure  14  is  a  circular  equatorial  orbit  with  a  semimajor 
axis  of  6678  km.  The  models  used  were  the  Harris -Priester,  Jacchia- Roberts,  Jacchia 
1960,  and  Jacchia  1971  models.  In  these  models,  variations  in  the  density  values  are  a 
result  of  the  different  density  data  collected  and  used  to  create  the  models.  AH  of  them 
show  the  basic  trends  in  density  with  time,  however  some  include  more  complex  inputs 
such  as  the  solar  flux  level  and  the  geomagnetic  activity.  The  geomagnetic  activity  refers 
to  the  strength  of  the  Earth’s  magnetic  field. 

The  diurnal  variations  are  clearly  defined  by  aU  four  models.  The  peaks  in  Figure 
14  were  found  to  correspond  to  approximately  1400  local  solar  time.  The  actual  values 
of  density  differ  between  the  models;  however,  they  are  on  the  correct  order  of  magnitude 


for  the  sateUite’s  altitude. 


Densltv 
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Time  (min) 


Figure  14:  Calculated  Densities  as  a  Function  of  Time  using  Four  Different 

Atmospheric  Models 
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Aiiiude  (Km) 

Figure  15:  Calculated  Densities  as  a  Function  of  Altitude  using  Four  Different 

Atmospheric  Models 

Figure  15  plots  density  vs.  altitude  for  the  same  case  in  Figure  14.  Figure  15 
demonstrates  the  similarity  among  the  four  models  where  the  increases  are  similar.  The 
Jacchia  1971  model  and  the  Jacchia-Roberts  model  closely  match  each  other  in  this  mn. 
The  similarity  is  because  the  Jacchia-Roberts  model  is  a  modification  of  the  Jacchia  1971 
model. 

Many  atmospheric  models  account  for  variations  in  the  solar  flux.  Figure  16  plots 
three  STK®  mns  using  the  Harris -Priester  model.  The  mns  differed  by  the  solar  flux 


input  into  the  propagator. 
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Figure  16:  Solar  Flux  Effects  on  Atmospheric  Density 

The  three  values  for  the  low,  average,  and  high  solar  flux  inputs  were  65  x  10“^^, 
150  X  10“  ,  and  250  x  10“  W/m  Hz  respectively.  Figure  16  shows  that  as  the  solar  flux 
increases,  the  density  at  a  given  altitude  increases  as  well. 

The  last  major  variation  in  density  taken  into  account  by  some  models  is  the 
effects  of  the  geomagnetic  activity.  A  higher  geomagnetic  index  causes  the  density  at  a 
given  altitude  to  increase,  while  a  lower  geomagnetic  activity  relates  to  a  lower  density  at 
a  given  altitude.  Figure  17  illustrates  three  mns  with  different  values  for  the  geomagnetic 
index.  The  geomagnetic  index  is  a  value  used  to  describe  the  activity,  where  a  higher 
number  indicates  the  geomagnetic  activity. 
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)(1o  '  Density  VS- Altitude 


Figure  17:  Geomagnetic  Activity  Effects  on  Atmospheric  Density 

The  density  calculation  algorithm  that  was  developed  demonstrates  the  complex 
variations  in  density  among  different  models.  Since  these  models  are  complex  and 
determining  density  values  involves  interpolations  among  different  tables,  it  is  difficult  to 
accurately  determine  the  error  in  the  calculated  density  and  reference  density.  However, 
based  on  the  errors  in  the  density  as  compared  with  the  1976  Standard  Atmosphere,  one 
can  expect  the  errors  to  increase  with  time. 
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Errors 

In  the  density  analysis,  there  are  several  sources  of  error.  These  errors  are  a  result 
of  uncertainties  in  equation  (9).  In  the  STK®  simulation,  values  of  position,  and  velocity 
are  calculated  as  a  result  of  the  propagation.  In  addition,  the  coefficient  of  drag  (Cd) 
remains  a  constant  value.  However,  reahstic  conditions  add  errors  into  these  values.  For 
example,  the  position  and  velocity  data  obtained  from  the  GPS  receiver  will  deviate  from 
the  hue  positions  resulting  from  errors  in  the  receiver’s  accuracy.  For  the  Ashtech  G- 12 
HDMA  GPS  receiver  board,  errors  of  3.0  m  CEP  and  6.0  m  in  altitude  can  be  expected. 

The  error  in  the  time  rate  of  change  of  the  semimajor  axis  is  low  compared  to 
other  possible  sources  of  error.  In  order  to  calculate  the  time  rate  of  change  of  the 
semimajor  axis,  I  used  MATLAB®  to  fit  a  polynomial  to  the  semimajor  axis  vs.  time  plot. 
A  second  order  polynomial  was  used  to  fit  five  data  point  sections  on  the  curve  of 
semimajor  axis  vs.  time.  By  taking  the  derivative  of  the  polynomial,  the  time  rate  of 
change  of  the  semimajor  axis  was  determined. 

When  random  errors  are  introduced  to  the  STK®  simulated  data,  errors  in  the 
calculated  density  are  greater  than  1%  and  are  deemed  unacceptable.  The  errors  are 
present  because  the  orbital  elements  calculated  from  the  position  and  velocity  data  vary 
greatly.  As  a  result,  the  derived  semimajor  axis  does  not  vary  consistently  with  a  satellite 
decaying  in  the  atmosphere.  Thus,  the  time  rate  of  change  of  the  semimajor  axis  due  to 
atmospheric  drag  is  not  accurately  represented.  When  processing  raw  data  from  a  GPS 
receiver,  a  smoothing  algorithm  is  necessary  to  eliminate  the  random  noise  in  the  GPS 
receiver.  Barring  effects  of  the  random  noise  on  the  semimajor  axis,  the  maximum 
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predicted  error  of  6  meters  in  altitude  will  result  in  an  error  in  density  by  approximately 

0.02%. 

Another  possible  source  of  error  is  the  coefficient  of  drag.  In  low- Earth  orbit,  the 
coefficient  of  drag  for  a  sphere  can  be  approximated  to  be  2.2,  however,  errors  in  this 
value  can  be  as  much  as  10%^^’^^.  The  variations  in  the  coefficient  of  drag  are  largely  a 
function  of  the  atmospheric  composition.  Thus,  a  10%  error  in  the  correct  value  of  the 
coefficient  of  drag  results  in  a  10%  error  in  the  calculated  density. 

Finally,  the  last  major  source  of  error  hes  in  the  motion  of  the  atmosphere.  The 
atmosphere  is  assumed  to  rotate  with  the  Earth,  with  rotational  velocity  bounds  of  zero 
and  the  Earth’s  angular  velocity.  At  lower  altitudes,  the  atmospheric  velocity  is  higher 
due  to  friction  and  as  the  altitude  increases,  the  velocity  decreases  towards  zero.  In  the 
1976  U.S.  Standard  Atmosphere  model,  the  entire  atmosphere  is  modeled  to  rotate  with 
the  Earth.  In  initial  mns  of  my  density  calculations,  I  assumed  that  the  velocity  of  the 
atmosphere  was  zero.  With  this  assumption,  the  maximum  errors  were  approximately 
12%.  Therefore,  errors  in  the  density  calculations  due  to  the  variations  in  the 
atmospheric  velocity  can  be  assumed  to  be  within  12%. 
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Spacecraft  Design 

The  driving  requirements  in  designing  a  spacecraft  to  carry  out  the  mission  of 
collecting  data  to  determine  atmospheric  density  are  the  payload  as  weU  as  the  shape  and 
mass.  The  primary  payload  is  an  accurate  GPS  receiver  to  determine  the  position  and 
velocity  of  the  sateUite.  However,  another  major  requirement  for  the  sateUite  is  to 
maintain  a  constant  baUistic  coefficient.  To  accomphsh  this,  the  sateUite’ s  mass  must  not 
change  and  its  cross  sectional  area  and  coefficient  of  drag  must  remain  constant  as  weU. 
Therefore,  the  best  configuration  for  the  spacecraft  would  be  a  spherical  shape  with  no 
propulsion  system.  This  would  ensure  that  the  baUistic  coefficient  remains  constant  in 
any  orientation. 

A  proposed  GPS  receiver  for  the  spacecraft  is  the  Ashtech  G- 12  High  Dynamics 
and  Missile  Applications  (HDMA)  receiver  board.  The  Ashtech  receiver  is  quaUfied  to 
orbital  velocities  as  weU  as  altitudes.  The  receiver  is  also  smaU  and  Ughtweight  making  it 
ideal  for  a  smaU  sateUite  design.  In  addition,  the  power  requirement  of  1.8  W  is  suitable 
for  sateUite  operations.  Estimated  accuracies  at  orbital  velocities  for  the  receiver  are  3.0 
m  CEP  for  position  and  0.05  m/s  for  velocity. 

The  mission  of  the  spacecraft  lends  itself  to  a  “store  and  forward” 
communications  architecture.  As  the  satelUte  orbits  the  Earth,  it  wiU  take  position  and 
velocity  readings  from  the  GPS  and  store  them  into  a  storage  device  on  the  spacecraft. 
Using  STK®,  readings  every  60  seconds  is  sufficient  to  determine  the  atmospheric 
density  in  the  orbit.  However,  readings  every  30  seconds  or  less  would  aUow  increased 
accuracy  and  filtering  of  random  noise.  Once  the  satelUte  passes  over  the  U.S.  Naval 
Academy  and  downloads  its  stored  data,  the  stored  position  and  velocity  can  be  erased  on 


the  spacecraft.  Once  downloaded  to  the  ground  station,  the  data  wiU  be  processed  to 
determine  the  densities  in  the  satellite’s  orbit. 
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Currently,  the  U.S  Naval  Academy  Small  Satelhte  program  is  in  the  process  of 
designing  the  Prototype  Communications  SateUite  or  PC- Sat.  This  sateUite  will  be 
testing  a  GPS  receiver  developed  by  the  German  Space  Operations  Center.  With  this 
payload,  position  and  velocity  data  can  be  retrieved  and  saved.  Although  the  ballistic 
coefficient  will  not  be  constant  due  to  its  shape,  the  opportunity  exists  for  testing  data 
processing  techniques.  Upon  launch  of  the  sateUite  designed  to  determine  atmospheric 
density,  data  processing  methods  can  be  in  place  to  calculate  atmospheric  density 
immediately. 
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Conclusions 

My  work  has  demonstrated  that  it  is  possible  determine  the  atmospheric  density  in 
a  satelhte’s  orbit  using  a  GPS  receiver.  Several  assumptions  are  necessary  to  yield 
accurate  results.  Most  of  these  assumptions  are  reasonable;  however,  the  assumption  that 
the  atmosphere  rotates  with  the  Earth  has  the  potential  for  creating  the  highest  error.  The 
atmosphere  has  a  velocity  that  can  range  from  zero  to  the  Earth’s  rotational  speed.  With 
an  accurate  model  of  the  atmospheric  winds,  it  would  be  possible  to  determine  density 
more  accurately.  Based  on  the  bounds  placed  on  the  atmospheric  winds,  the  greatest 
possible  error  in  my  algorithm  is  12%.  This  error  is  too  large  to  allow  reasonable 
calculations  of  density  and  therefore,  future  modifications  must  be  made. 

Eor  my  algorithm  to  be  suitable  for  the  USNA  small  sateUite  program,  my 
algorithm  should  include  a  filter  to  eliminate  the  noise  from  the  GPS  data.  Even  small 
inaccuracies  in  the  GPS  data  cause  large  errors  in  the  density  calculations.  Another 
improvement  to  my  algorithm  would  be  to  eliminate  the  need  for  SateUite  Tool  Kit®.  By 
programming  a  propagator  into  the  algorithm,  it  wiU  simpUfy  the  density  analysis 
because  frequent  updates  of  the  orbital  elements  could  be  done  automaticaUy.  In 
addition,  the  orbital  element  updates  could  be  done  more  frequently,  driving  the  errors 
down  significantly.  This  project  has  laid  the  groundwork  for  the  next  generation  of  the 
USNA  smaU  satelUtes.  With  the  suggested  modifications,  fiiture  Midshipmen  can  use  the 
algorithm  to  create  a  database  for  scientists  to  use  in  creating  a  new,  more  accurate 


atmospheric  density  models. 
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Appendix  A:  Glossary 

Argument  of  Periapsis:  The  angle  measured  between  the  hne  of  nodes  and  the  eccentricity 
vector. 

Ballistic  Coefficient:  A  term  used  in  aerodynamics  defined  as  the  mass  divided  by  the  product 
of  the  coefficient  of  drag  and  the  cross  sectional  area 

Circular  Error  Probable:  The  radius  of  a  circle  where  50%  of  the  computed  positions  will  he 
within  the  circle. 

Circular  Orbit:  An  orbit  where  the  eccentricity  of  the  orbit  is  equal  to  zero. 

Coefficient  of  Drag:  A  dimensionless  value  which  reflects  the  sateUite’s  susceptibihty  to  drag, 
forces.^  (Cd) 

Direct  Orbit:  An  orbit  with  an  inchnation  less  than  90°. 

Earth’s  Gravitational  Parameter  The  term  that  is  defined  by  the  product  of  the  gravitational 
constant  and  the  Earth’s  mass  (|a.=  G'*m^). 

Eccentricity:  The  parameter  of  an  eUipse  defining  the  shape  of  the  orbit,  which  is  the  ratio  of  the 
distance  between  the  two  foci  to  the  major  axis  (2a). 

Eccentricity  vector:  A  vector  that  points  in  the  direction  of  periapsis  and  has  the  magnitude  of 
the  eccentricity. 

Eccentric  Orbit:  An  orbit  with  an  eccentricity  between  0  and  1.  An  eccentricity  of  0  defines  a 
circular  orbit  while  an  eccentricity  of  1  defines  a  parabohc  orbit. 

Flight  path  angle:  The  angle  between  the  velocity  vector  and  the  local  horizon,  where  the  local 
horizon  is  perpendicular  to  the  position  vector. 

Inclination:  The  angle  between  an  orbit’s  angular  momentum  vector  and  the  Earth’s  angular 
momentum  vector.  It  can  also  be  measured  by  the  angle  between  the  orbital  plane  and  the 
Earth’s  equatorial  plane. 

Line  of  nodes:  The  vector  connecting  the  center  of  the  Earth  to  the  ascending  node.  The 
ascending  node  is  the  point  where  the  satelhte  crosses  the  equatorial  plane. 

Longitude  of  Ascending  node :  The  angle  between  the  /  unit  vector  (in  the  Geocentric 
Equatorial  coordinate  system)  and  the  hne  of  nodes,  where  the  /  vector  points  in  the  direction  of 

‘  Vallado,  498. 
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the  first  point  of  Aries.  The  first  point  of  Aries  is  defined  by  the  hne  from  the  Earth  to  the  Sun, 
when  the  Earth  is  at  the  vernal  equinox. 

Low-Earth  Orbit:  An  orbit  with  altitudes  below  800  km. 

Mean  Motion:  The  mean  angular  rate  of  a  sateUite’s  orbit  defined  by  «  = 

Periapsis:  The  point  of  the  satellite’s  orbit  that  is  the  closest  to  the  central  body. 

Perigee:  The  point  of  the  sateUite’s  orbit  that  is  the  closest  to  the  Earth. 

Polar  Orbit:  An  orbit  with  an  incUnation  of  90°. 

Retrograde  orbit:  An  orbit  with  inclinations  greater  than  90°. 

Scale  Height:  The  distance  where  the  atmospheric  density  decreases  by  a  factor  of  - . 

e 

Semimajor  Axis:  Half  of  the  distance  of  the  eUipse’s  length. 


True  Anomaly:  The  angle  between  the  eccentricity  vector  and  the  position  vector  as  shown  in 
Eigure  1  on  page  11. 


Appendix  B:  Additional  Results  -  Density  /  Error  Plots 
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Figure  B.3  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.5  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.7  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.9  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.ll  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.13  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.15  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.17  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
a  =  6678  km 
e  =  0.001 
i  =  0° 


Density  (kg/tn  ) 


B-20 


Figure  B.19  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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e  =  0.001 
i  =  45° 


Density  (kg/tn  ) 


B-22 


Density  vs.  Altitude 


Figure  B.21  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.23  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.25  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
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Figure  B.27  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
a  =  6778  km 
e  =  0 
i  =  45° 
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Figure  B.29  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
a  =  6778  km 
e  =  0.001 
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Figure  B.31  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
a  =  6778  km 
e  =  0.001 
i  =  45° 
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Figure  B.33  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
a  =  6778  km 
e  =  0.01 
i  =  0° 
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Figure  B.35  Density  vs.  Altitude  For  an  Orbit  Defined  By: 
a  =  6778  km 
e  =  0.01 
i  =  45° 


Appendix  D:  Gaussian  Variation  of  Parameters 
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Vallado*  derives  equation  (7)  from  the  Gaussian  form  of  Lagrange’s  Variation  of 
Parameters.  The  Gaussian  form  defines  the  changes  in  the  classical  orbital  elements  as  a 
function  of  disturbing  forces.  The  equations  are  defined  in  terms  of  Fr  and  Fs,  the  components 
of  a  force  F.  These  components  are  in  the  RSW  coordinate  system.  The  R  axis  points  in  the 
radial  direction.  The  S  axis  is  in  the  orbital  plane  and  is  perpendicular  to  the  R  axis,  pointing  in 
the  direction  of  motion.  The  W  axis  is  perpendicular  to  both  the  R  and  S  axes,  creating  a  right 
hand  coordinate  system. 


da  2e  *  sin^ ) 


(D.l) 


The  flight  path  angle  can  be  used  to  divide  the  force  due  to  drag  into  components  in  the  RSW 
coordinate  system  as  seen  in  equations  (D.2)  and  (D.3). 


(D.2) 


1  C  *  A 

Fs^--P  — - *COS((t)  fpa) 

2  m 


(D.3) 


The  flight  path  angle  is  defined  as  the  angle  between  the  velocity  vector  and  the  local 
horizon.  The  sine  and  cosine  of  the  flight  path  angle  can  be  written  as  a  function  of  the 
eccentricity  and  true  anomaly  as  seen  in  equations  (D.4)  and  (D.5). 


sin((l)/pj  =  ■ 


e*sin(v) 
ll  +  l* e*  costs' ) 


(D.4) 


cos((l)ypJ  =■ 


1+  e*cos(v) 

/l  +  2  *  e*  cos(v) 


(D.5) 


Substituting  equations  (D.2)  and  (D.3)  into  equation  (D.l),  VaUado  obtains: 


da  2e*smty)(  1  ^ 


nvl  -  e 


Vrel  *sin((l)^^J+- 


f  1  Cn*A  2  ) 

=  --9— - ^rel  *COS(. 

2  2  m 


(D.6) 


By  reducing  equation  (D.6),  Vallado  obtains  equation  (D.7) 


Cn*A 


(e*sin(v))^ 


p{l  +  e*  cos(v)) 


nyl  —  e^  \  -Jl  + +2*e*costy)  rJl  + +2* e*costy) 


(D.7) 


Further  simplification  yields: 


D-2 


da  A  2 

( 

;  1 
l-l-e  -l-2*e*cos(v') 

(e*sin^))^  -l-(l-l- e*cos(v'))^ 

j  P  ^rel 

at  m 

1 

rfsll  —  e^ 

) 

1-1-  -1-  2*e*cos(v') 

Finally,  equation  (D.8)  can  be  reduced  to  obtain: 


da  C[)*A  2 

-r=“P - 

dt  m 


i 


n'Jl  —  e" 


(D.8) 


(D.9) 


*  David  A.  Vallado,  Fundamentals  of  Astrodvnamics  and  Applications.  New  York:  The  McGraw  Hill  Companies, 
1997, 559-568, 604-605. 
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Appendix  E:  Modeling  Atmospheric  Rotation 

King-Hele^  derives  a  scalar  equation  for  calculating  the  relative  velocity  of  a  satellite 
with  respect  to  the  atmosphere.  The  assumption  made  in  this  derivation  is  that  the  atmosphere 
rotates  with  the  Earth. 

King-Hele  begins  with  the  vector  form  for  the  relative  velocity  of  the  sateUite  with 
respect  to  the  atmosphere  where  is  the  relative  velocity,  v  is  the  inertial  velocity  of  the 

sateUite,  and  is  the  inertial  velocity  of  the  atmosphere: 

^rel=^-^A  (E-1) 

By  using  the  law  of  cosines,  the  scalar  values  can  be  used  to  find  where  y  is  the  angle 
between  v  and  . 


vli^v^  +  V^-2*v*V^*cos(/)  (E.2) 

Assuming  the  atmosphere  rotates  with  the  Earth,  can  be  written  as  equation  (E.3)  where  r  is 
the  magnitude  of  the  sateUite’s  position  vector,  is  the  angular  velocity  of  the  Earth,  and 
is  the  latitude  of  the  sateUite. 

VA  =  r*^(B*^o^(^SAT)  (E-3) 

Substituting  equation  (E.3)  into  equation  (E.2)  yields: 

=v^  +  r^  *cos((;|>j^y.)*cos(Y)  (E.4) 

Using  relationships  between  a  spherical  triangle  formed  by  the  incUnation,  latitude  of  the 
sateUite,  and  the  angle  between  the  velocity  vectors  results  in  the  foUowing: 

COS(0  =  COS((|>j^^)  *cos(y)  (E.5) 


The  resulting  equation  in  terms  of  the  incUnation  and  the  satelUte’s  latitude  can  be  written  as: 


2  2 
Vrel  =  V 


1-- 


-COS(0 


+  r^*Q.(D^  *{cos^{<\>SAT)-cos^ii)) 


(E.6) 


*  Desmond  King-Hele,  Satellite  Orbits  in  an  Atmosphere:  Theory  and  Applications.  London,  England:  Blackie  and 
Son,  Ltd,  1987,  29-30. 
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